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Several extensions of General Relativity and high energy physics include scalar fields as extra 
degrees of freedom. In the search for predictions in the non-linear regime of cosmological evolution, 
the community makes use of numerical simulations in which the quasi-static limit is assumed when 
solving the equation of motion of the scalar field. In this Letter, we propose for the first time a 
method to solve the full equations of motion for scalar degrees of freedom. We run cosmological 
simulations which track the full time and space evolution of the scalar field, and find striking 
differences with respect to the commonly used quasi- static approximation. This novel procedure 
reveals new physical properties of the scalar field and reveals concealed astrophysical phenomena 
which were hidden in the old approach. 



General Relativity (GR) is the foundation stone of the 
standard cosmological model. The assumption that GR 
describes gravity at cosmological scales leads to two ex- 
tra building blocks of the model: dark matter and dark 
energy. In spite of the model's ability to match a large 
number of observations, the nature of the two dark com- 
ponents is still unclear. Several extensions of General 
Relativity and the Standard Model of Particle Physics 
have been proposed to explain the dark sector. Most of 
this extensions include scalar fields as extra degrees of 
freedom which interact with matter and lead to modifi- 
cations to standard gravity As GR is proven to be 
valid in solar system scales, any modification introduced 
must fulfil the requirement of reducing to Einstein grav- 
ity at small scales. This is assured through screening 
mechanisms [2|-l4j. Such processes are highly non-linear, 
therefore are better probed in the non-linear regime of 
cosmological structure formation in the scales of galaxies 
and clusters of galaxies. 

With the intention to probe the nature of the dark sec- 
tor, modifications to GR and the indispensable screening 
mechanisms, the community makes use of N-body sim- 
ulations. The main assumption behind such simulations 
is the quasi-static limit when solving the equation of mo- 
tion of the extra scalar degrees of freedom (i.e. neglect 
time derivatives) [5-8]. The validity of this quasi-static 
hypothesis was never tested. 

In this Letter, we propose for the first time a method 
to solve the full equations of motion for scalar degrees of 
freedom. We run cosmological simulations which track 
the full time and space evolution of the scalar field, and 
find striking differences with respect to the commonly 
used quasi-static approximation. Our novel procedure 
reveals new physical properties of the scalar field and 
reveals concealed astrophysical phenomena which were 
hidden in the old approach. We show how such results 
may have an impact on both the theoretical studies of 
Modified Gravity and its astrophysical probes. 

We focus on a general class of scalar-tensor theories 
which present the symmetron screening mechanism 



However, our results and techniques can be straightfor- 
wardly generalised to any theory with a scalar degree of 
freedom, independently of the screening mechanism that 
is invoked. 

The symmetron model is defined by the potential: 
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which leads to the following equation of motion: 
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where \i and M are mass scales, p is the matter density. 
We normalize the field <p with its vacuum expectation 
value 0o = ^- . Defining x = ^ > an d taking into account 
that the background density acquires a value Pssb — 
M 2 p? at zssb, we obtain: 
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where rj is the matter density field normalised with the 
background density, and Ao = ^7=^ is the range for the 
field that corresponds to 77 = 0. 

The Newtonian potential 0at is given by: 
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where 5 is the over-density denned as Sp/pb, Pb is the 
background density, 0at is the perturbation in the metric. 

The dark matter component can be tracked by the par- 
ticles that follow the geodesies equation: 
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The standard way of solving equation (j3j) is to assume the 
quasi-static limit and use a standard multigrid method to 
solve the resulting elliptic equation. In this Letter we go 
beyond this approach and solve the complete equation. 
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Our new method is based on the fact that the equation 
of motion for the scalar field, eq.(j3]), is formally equiva- 
lent to the geodesies equation, eq.flSJ). We can, therefore, 
apply a leap-frog scheme, not to positions and velocities 
of particles as in an N-body code, but to the scalar field 
X and its time derivative on a grid. The change q = a 3 x 
gives the following first order equations: 



Q 
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By using second order discretization in time and imple- 
menting a leap frog scheme, we obtain the following equa- 
tions for the field and the particles for the time step n: 
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where is the total potential that includes GR and the 
scalar field and a is employed as time variable. This 
method was successfully applied to compute solutions of 
the growth equation for linear density perturbations cou- 
pled to a non-linear Poisson equation [9]. Furthermore, 
the coupling with the geodesies was employed to generate 
initial conditions for non-linear gravities Q. 

The implementation of the solution requires initial con- 
ditions for the variables x an d q. We choose to fix x — 
at the initial time, which corresponds to a fully screened 
scalar field as is expected at early times. For g, we use a 
small random number. 

Our new method and its parallelisation were imple- 
mented in the N-body code presented in [9] , to which the 
leap-frog associated to the time evolution of the scalar 
field was added. The algorithm and the 3D solver were 
tested satisfactorily against ID analytic linear solutions 
that were obtained for a fixed density distribution given 
by a gaussian profile located in the center of the box. 

For the first time in cosmological simulations of modi- 
fied gravity one can compute the speed at which the new 
scalar degree of freedom travels. To this end, we use 
the most straightforward definition, which is the velocity 
at which perturbations propagate in space away from a 
given source. The origin of the perturbations is given by 
a fixed ID Gaussian density profile with a dispersion of 
1 Mpc/h located in the center of the domain. The ini- 
tial conditions are x — 1 an d q — 0, which means that 
we start the simulations with an unperturbed vacuum 



100 



Q. o 



-100 




0.9 



0.8 



0.7 



-150 



0.92 0.925 0.93 0.935 0.94 0.945 0.95 0.955 0.96 



FIG. 1: Normalised scalar field % as a function of space and 
time for a fixed Gaussian density distribution. For clarity, the 
unperturbed regions where % = 1 are shown in black. The 
green lines are null geodesies that pass through the center of 
the halo at the initial time. Here, zssb — 1 and Ao = 1 
Mpc/h. 



field and measure the speed of the perturbations induced 
by the density as they move away from the source. A 
graphical representation of the speed of sound can be 
seen in FigJTJ where we show x as a function of space 
and time. The bulk of the mass that induce the oscilla- 
tions is located in the center of the domain. The green 
lines are null geodesies that depart from the center of the 
halo at the initial time. The figure shows that the front 
wave of the perturbations in the scalar field moves with 
the same speed as light, as expected due to its standard 
kinetic term. We performed several simulations with dif- 
ferent configurations for the underlying density distribu- 
tions and found no deviation from this value for ranges 
A from 1(T 4 to 1 Mpc/h. 

In order to highlight the consequences of taking into 
account non-static terms in the cosmological evolution of 
the scalar field, we present results obtained from a cos- 
mological simulation that was run with a box at redshift 
z = of 128 Mpc/h and using 128 3 particles. The scalar 
field and gravitational potential were tracked in a uniform 
grid with 128 nodes per dimension. As the evolution of 
the matter density and the scalar field have different time 
scales and the aim of this Letter is to understand only the 
evolution of the scalar field, we evolve the matter com- 
ponent by taking into account the GR equations. The 
small differences that will appear in the matter compo- 
nent because of considering also the fifth force in its evo- 
lution are not expected to change the main properties of 
the non- static solution that we intend to study here. The 
background evolution is given by a flat ACDM cosmology 
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FIG. 2: The four color panels show the normalised scalar field x f° r a particular slice of the simulation at four different redshifts 
(z = 0.807, 0.366, 0.03 and 0.008). The white lines are density contours. Up- left: scalar field developing in voids (zssb = 0.5). 
Up-Right: a domain wall which is stable from its formation up to z ~ 0.03. Bottom-left: collapse of the wall. Bottom-right: a 
spherical wave originated in that collapse and a dark blob (see the upper- left region where the wall was). The extra panel to 
the right shows a cut of the last color panel through a row with y = 108 Mpc/h (the curves are density and normalised scalar 
field). Note the dark blob locate in x ~ 25 Mpc/h. 



(Q m = 0.267, tt A = 0.733 and H = 71.9 km/sec/Mpc). 
The initial conditions for matter were generated using 
Zeldovich approximation with standard gravity with the 
package Cosmics [10j . The particular symmetron param- 
eters employed are zssb = 0.5 and Ao = 1 Mpc/h, which 
provide a cosmological evolution close to ACDM [H| . 

The time steeping is uniform in a space during the 
simulations. We used 10 4 time steps to evolve GR and 
the position of the particles. The evolution of the scalar 
field has a much shorter time scale than GR and thus, 
it must be tracked with a much shorter time step. We 
made runs with 9.6 x 10 5 and 1.92 x 10 6 time steps and 
found the same statistical distribution of the scalar field 
(shown in Figj3]), which shows that this number of steps 
is enough for the spatial resolution employed. 

Fig|2] shows the distribution of the non-static scalar 
field for a given slice of the box at four different redshifts 
(z = 0.807,0.366,0.03 and 0.008). The white lines are 
contours of the matter density. The different snapshots 
were chosen to highlight different properties of the non- 
static solutions. 

The symmetron has an effective potential which ac- 



quires two minima when the symmetry is broken and 
thus, the scalar field can adopt different signs in different 
regions in space. In the static cosmological simulations 
one needs to fix the sign of the field to be unique all over 
the box (relaxing this constraint make iterative solvers to 
fail to converge). As in our simulations the evolution of 
the scalar field is followed in a self consistent way, we can 
study the validity of this constraint. We find that the 
sign of x i s indeed not necessarily unique. On the con- 
trary, after zssb, the box appears divided in different 
domains in which the scalar field chooses different signs. 
We find that the domain walls associated to this decom- 
position can be stable during long period of time until 
they collapse. The upper panels of the Fig|2] show the 
formation of a domain wall in the upper-left part of the 
box. This particular wall does not survive until redshift 
z = 0. In the bottom- left panel, it is possible to see the 
wall collapsing. The energy that was contained in it is re- 
leased in the form of scalar waves that can travel several 
megaparsecs. The front wave associated to this particu- 
lar event can be found in the lower-right panel crossing 
the center of the box. It is important to note that in the 
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FIG. 3: Histograms of normalised scalar field \ on the grid 
for three different redshifts of the cosmological simulation. 



static case, the matter distribution only fixes the bound- 
ary conditions. However, the exact time evolution of \ 
is a function not only of the boundary conditions, but 
also of its history, which is not described by the static 
equation. Even if the domain walls can be obtained from 
the static solution, the static equation does not provide 
enough information to describe their evolution in time, 
neither the release of energy and propagation of scalar 
waves. All these factors affect the behaviour (range and 
strength) of the fifth force induced by the scalar field and 
so the formation of cosmological structures. 

We run several simulations with different initial seeds 
for the realisations of the initial density distribution. We 
find that the formation of domains is a common phe- 
nomenon, as it is that not all of them survive until red- 
shift z = 0. Therefore, observational signatures of this 
process are more likely to be found at redshifts closer to 
symmetry breaking. The question whether our universe 
has a unique sign can be only answered by running larger 
simulations. 

The distribution of signs for the simulation presented 
in this Letter can be seen in Figj3j where we show his- 
tograms of the scalar field made over the whole grid at 
different redshifts. At redshift close to symmetry break- 
ing, the scalar field is almost fully screened and thus, it 
has a large peak in zero. Later on, two large picks de- 
velop in the distribution, which corresponds to domains 
with positive and negative \. For the particular simu- 
lation shown here, the bimodal distribution subsist until 
z = 0. We find that this is not a general result, but 
that by changing the initial seed of the initial conditions 
it is possible to obtain distributions completely skewed 
towards positive or negative values. 



The non-linearity of modified gravity can induce a de- 
coupling between dynamical and real mass, i.e. the mod- 
ified force can converge towards regions devoid of matter 
. Here we study the possibility that non-static terms 
can also induce similar effects and find that the answer is 
indeed positive. We invite the reader to focus his/her at- 
tention in the lower-right color panel of Fig|2j The figure 
shows that each feature in the scalar field distribution 
has a counterpart in the density distribution with excep- 
tion of the minimum of |x| at (x,y) ~ (25,108) Mpc/h 
(it can be seen as a black region in the color figure). We 
find that this feature of the solution of the field equation 
is not a consequence of an external field as in , but a 
consequence of the non-static terms: the effect does not 
appear in the static solution that correspond to the same 
density distribution. Furthermore, the depression in x 
moves in space, a characteristic that cannot be explained 
by using the static approximation. We find the same ef- 
fect in several simulations, usually associated to domain 
walls. While the energy associated with the scalar field is 
not enough to explain the dark blobs that are present in 
several clusters of galaxies 12-14|, this peculiar charac- 
teristic of the non-static solution must be taken into ac- 
count when developing new observables to detect scalar 
fields and probe modifications to GR. 

In summary, we present a new method to solve the full 
equation of motion for scalar-tensor theories that does 
not rely on the quasi-static approximation. The method 
is based on a leap-frog scheme and can be generalised 
in a trivial way to a number of commonly studied mod- 
ified gravity models by changing the source function in 
eq.0. In order to test the consequence of the non-static 
terms on the evolution of the scalar field, we run cosmo- 
logical simulations with our new solver. We find several 
new striking physical features which are not present in 
the usual static approach. For instance, the scalar de- 
gree of freedom develops domain walls whose evolution 
in time cannot be described by the static equation. The 
collapse of these domain walls are source of scalar waves 
that travel across several Mpc and may leave a cosmo- 
logical imprint in galaxies and clusters at those scales. 
Finally, our algorithm shows that it is possible to decou- 
ple the potential wells in the scalar field from overdense 
regions, leading to the decoupling of real and dynamical 
mass. The questions related to the impact that these 
large differences that exist between static and non-static 
solutions has in the matter distribution itself is left for 
future work. Simulations show that matter can fall in 
the domain walls producing a local increase in the mat- 
ter density. This makes the domain walls to leave an 
imprint in the cosmic shear that should be possible to 
observe with large lensing surveys such as the incoming 
ESA Euclid survey. Furthermore, as this modifications 
in the density depends of the position of the walls and 
not of the matter distribution, the non-static terms pro- 
vide a new source of non-gaussianities that should be 
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taken into account. While the predictions related to do- 
main walls are specific of the symmetron model, there are 
other observables related to the oscillations in the scalar 
field that are model independent. In this sense, studies of 
mass accretion in galaxies for instance through L a emis- 
sion (which occur in regions where the scalar fields oscil- 
late) and generation of X-rays in shock waves induced by 
the oscillations will be crucial when searching for tighter 
constraints for scalar-tensor theories. The implications 
of these effects if observed in near future astrophysical 
probes at galactic and cluster scales would open a new 
window to test novel gravity theories and high energy 
models with cosmological simulations of the nonlinear 
regime. 



[1] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, 

Phys.Rep. 513, 1 (2012), 1106.2476. 
[2] K. Hinterbichler and J. Khoury, Physical Review Letters 

104, 231301 (2010), 1001.4525. 
[3] J. Khoury and A. Weltman, Phys. Rev. Lett. 93, 171104 



(2004), astro-ph/0309300. 
[4] T. S. Koivisto, D. F. Mot a, and M. Zumalacarregui, 

Physical Review Letters 109, 241102 (2012), 1205.3167. 
[5] A.-C. Davis, B. Li, D. F. Mota, and H. A. Whither, ArXiv 

e-prints (arXiv:1108.3081) (2011), 1108.3081. 
[6] B. Li, G.-B. Zhao, R. Teyssier, and K. Koyama, JCAP 

1201, 051 (2012), 1110.1379. 
[7] M. Baldi, Physics of the Dark Universe 1, 162 (2012). 
[8] F. Schmidt, Phys.Rev. D80, 043001 (2009), 0905.0858. 
[9] C. Llinares, PhD thesis. ISBN: 978-90-367-4760-8 

http: / / dissertations. ub. rug. nl/faculties/science/201 1 / c. llinares 

(2011). 

[10] E. Bertschinger, ArXiv Astrophysics e-prints (1995), 

arXiv:astro-ph/9506070. 
[11] A. Knebe, C. Llinares, X. Wu, and H. S. Zhao, Astrophys. 

J. 703, 2285 (2009), 0908.3480. 
[12] M. J. Jee, A. Mahdavi, H. Hoekstra, A. Babul, J. J. 

Dalcanton, P. Carroll, and P. Capak, Astrophys. J. 747, 

96 (2012), 1202.6368. 
[13] A. Mahdavi, H. Hoekstra, A. Babul, D. D. Balam, and 

P. L. Capak, Astrophys. J. 668, 806 (2007), 0706.3048. 
[14] S. Riemer-Sorensen, K. Pedersen, S. H. Hansen, and 

H. Dahle, Phys. Rev. D 76, 043524 (2007), arXiviastro- 

ph/0610034. 



